Freshwater Mulit-Stressor Synthesis - Figures

This R notebook contains the code used to create the figures in the paper: “Interactions among anthropogenic stressors in freshwater ecosystems: a systematic review of 2,396 multiple-stressor experiments” by James A. Orr, Samuel J. Macaulay, Adriana Mordente, Benjamin Burgess, Dania Albini, Julia G. Hunn, Katherin Restrepo-Sulez, Ramesh Wilson, Anne Schechner, Aoife M. Robertson, Bethany Lee, Blake Stuparyk, Delezia Singh, Isobel O’Loughlin, Jeremy J. Piggott, Jiangqiu Zhu, Khuong V. Dinh, Louise C. Archer, Marcin Penk, Minh Thi Thuy Vu, Noël P.D. Juvigny-Kkenafou, Peiyu Zhang, Philip Sanders, Ralf B. Schäfer, Rolf Vinebrooke, Sabine Hilt, Thomas Reed, Michelle C. Jackson.

Set-up

• Set up R enviornment

# R version 4.3.1 (2023-06-16)

rm(list = ls())           # clear the environment 

# Packages
library(sf)               # classes and functions for vector data
library(raster)           # classes and functions for raster data
library(tidyverse)        # organising and manipulating data
library(spData)           # for world shapefile 
library(ggVennDiagram)    # for venn diagram figure
library(igraph)           # to create co-occurrence network 
library(knitr)            # to create nice tables
library(reshape2)         # for the melt function 
library(magick)           # for making gifs 
library(ggstream)         # for proportion plots
library(gridExtra)        # for arranging multiple ggplots

• Load data

The experiments dataset contains 2396 rows, one for each experiment. The combinations dataset contains 4712 rows, one for each combination of stressors, which we use for the co-occurrence analysis. The SpecificCombinations column in the experiments dataset contains the coding used to create the combinations. A description of all of the columns can be found in the supporting information of the paper.

Figure 1 - Overview

• Map of experiments

# all lower case first
experiments$Country <- tolower(experiments$Country)

# organise experiment data
countries <- experiments %>%
  mutate(counts = 1) %>%
  group_by(Country) %>%
  summarise(experiments = sum(counts)) %>%
  rename(name_long = Country)

# Mercator - epsg: 3857
# Mollweide - `st_crs('ESRI:54009')`
# Robinson - `st_crs('ESRI:54030')`
# Winkel Triple - `st_crs('ESRI:54019')`

# load country shapes
world <- st_transform(spData::world, st_crs('ESRI:54019')) 

# cleaning to make the dataframes consistent
world$name_long <- tolower(world$name_long)
countries$name_long <- tolower(countries$name_long)
countries$name_long <-replace(countries$name_long,
                           countries$name_long=="russia", "russian federation")
countries$name_long <-replace(countries$name_long,
                           countries$name_long=="the netherlands", "netherlands")
countries$name_long <-replace(countries$name_long,
                           countries$name_long=="uk", "united kingdom")
countries$name_long <-replace(countries$name_long,
                           countries$name_long=="usa", "united states")
countries$name_long <-replace(countries$name_long,
                           countries$name_long=="slovak republic", "slovakia")

# join experiment data to countries 
world <- world %>%
  left_join(countries, by = "name_long") %>%
  select(name_long, experiments)

• Venn diagram of responses

# List of items
x <- list(Physiological = which(experiments$Physiological == "yes"),
          Individual = which(experiments$Individual == "yes"), 
          Population = which(experiments$Population == "yes"),
          Community = which(experiments$Community == "yes"),
          Ecosystem = which(experiments$Ecosystem == "yes"))

# Venn diagram
ggVennDiagram(x, label = "count", label_alpha = 0, edge_size = 0,
              label_size = 3,
              category.names = c("Phy.",
                                 "Ind.",
                                 "Pop.", 
                                 "Com.",
                                 "Eco.")) +
  
  labs(fill = "experiments") +
  
  scale_fill_gradient(low = rgb(0.8, 0.8, 0.8), 
                      high = rgb(0.4, 0.4, 0.4), trans = 'sqrt') +
  
  scale_color_manual(values = c("grey30", "grey30", "grey30", "grey30", "grey30"))

  
rm(x)

• Number of Stressors and Fully Factorial

Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
Please use the `linewidth` argument instead.

• Stressor Levels

• System and Habitat

• Duration

Warning: NAs introduced by coercion

• Taxonomic Groups

• Temporal Aspects

• Publication growth rate

• Compared to “Biology” in WOS

Joining with `by = join_by(Year)`

# Define the year intervals
year_intervals <- c(1970, 1980, 1990, 2000, 2010, 2020)

# Calculate annual growth rates for experiments and WOS.Biology within each interval
annual_growth_rates <- data.frame()

for (i in 1:(length(year_intervals) - 1)) {
  start_year <- year_intervals[i]
  end_year <- year_intervals[i + 1]
  
  interval_growth_rates <- experiments_year %>%
    filter(Year >= start_year, Year < end_year) %>%
    summarise(experiments_growth = (((last(experiments)/first(experiments))^(1/10)) - 1) * 100,
              wos_growth = (((last(WOS.Biology)/first(WOS.Biology))^(1/10)) - 1) * 100)
  
  annual_growth_rates <- bind_rows(annual_growth_rates, interval_growth_rates)
}

# Extract and return the annual growth rates as vectors
experiments_growth_rates <- annual_growth_rates$experiments_growth
wos_growth_rates <- annual_growth_rates$wos_growth

experiments_growth_rates

[1]      Inf      Inf 25.24845 12.33498 10.95004

wos_growth_rates

[1] 5.702093 1.105205 3.182726 2.992679 2.682283

rm(i, start_year, end_year, year_intervals, all_years, 
   wos_growth_rates, interval_growth_rates, WOS)

Table 1 - Taxonomy

identity_data <- experiments %>%
  select(contains("Identity")) %>%
  gather(value = "Identity")

class_data <- experiments %>%
  select(contains("Class")) %>%
  gather(value = "Class")

taxonomy <- data.frame(Class = class_data$Class, Identity = identity_data$Identity)

# Filter out rows with NA in the "Class" column (e.g., StressorTenClass is mostly NA)
taxonomy <- taxonomy %>%
  filter(!is.na(Class))

# Group by the "Class" column and summarize the unique identities
taxonomy_table <- taxonomy %>%
  group_by(Class, Identity) %>%
  summarize(Count = n()) %>%
  ### option to select by frequency of abundances ###
  #filter(Count >= 5) %>%
  group_by(Class) %>%
  ### where the magic happens ###
  summarize(Unique_Identities = paste0(Identity, " (", Count, ")", collapse = ", "), .groups = "drop")

`summarise()` has grouped output by 'Class'. You can override using the `.groups` argument.

# Reorder the rows based on the desired order for the taxonomy table
desired_order <- c("habitat alteration", "hydrology",  
                   "radiation", "sound", "temperature",
                   "UV light", "visible light", "water clarity",
                   
                   "microplastics", "nanoparticles", 
                   
                   "acidity", "alkalinity", "antibiotic", "carbon dioxide", 
                   "fungicide", "general biocide", "herbicide", "insecticide", 
                   "metals", "nutrients", 
                   "other chemicals", "oxygen", "pharmaceuticals", "salinity",
                   
                   "cyanotoxin",
                   
                   "biological alterations", "biological control", "disease",
                   "domestic species", "non-native species",
                   
                   "composite stressors")

taxonomy_table <- taxonomy_table[match(desired_order, taxonomy_table$Class), ]

# Create a nicely formatted table
#kable(taxonomy_table, format = "markdown", col.names = c("Class", "Unique Identities"))

# Taxonomy dataset to use for analyses 
taxonomy <- taxonomy %>%
  group_by(Class) %>%
  summarize(N = n()) %>%
  mutate(Nature = ifelse(Class %in% c("habitat alteration", "hydrology",  
                   "radiation", "sound", "temperature",
                   "UV light", "visible light", "water clarity"),
                   "Physical", 
                   ifelse(Class %in% c("biological alterations", 
                                       "biological control", "disease",
                                       "domestic species", "non-native species"),
                   "Biological",
                   ifelse(Class == "composite stressors", 
                          "Mixed", 
                        ifelse(Class == "cyanotoxin", "Biological-Chemical",
                               ifelse(Class %in% c("microplastics", "nanoparticles"), 
                                      "Chemical-Physical",
                                      "Chemical"))))))

Figure 2 - Stressor Classes

• Stressor natures over time

class_data <- experiments %>%
  select(c(Year, contains("Class"))) %>%
  pivot_longer(cols = contains("Class")) %>%
  drop_na() %>%
  select(-name) %>%
  mutate(Nature = ifelse(value %in% c("habitat alteration", "hydrology",  
                   "radiation", "sound", "temperature",
                   "UV light", "visible light", "water clarity"),
                   "Physical", 
                   ifelse(value %in% c("biological alterations", 
                                       "biological control", "disease",
                                       "domestic species", "non-native species"),
                   "Biological",
                   ifelse(value == "composite stressors", 
                          "Mixed", 
                        ifelse(value == "cyanotoxin", "Biological-Chemical",
                               ifelse(value %in% c("microplastics", "nanoparticles"), 
                                      "Chemical-Physical",
                                      "Chemical"))))))

### merge all years pre 1991 into 1990 (so that kernel density estimation is more accurate)
class_data <- class_data %>%
  mutate(Year = ifelse(Year < 1991, 1990, Year))


### summarise by Nature 
nature_data <- class_data %>%
  group_by(Year, Nature) %>%
  summarise(Count = n())

`summarise()` has grouped output by 'Year'. You can override using the `.groups` argument.

### summarise by Class 
class_data <- class_data %>%
  group_by(Year, value) %>%
  summarise(Count = n())

`summarise()` has grouped output by 'Year'. You can override using the `.groups` argument.

nature_data$Nature <- factor(nature_data$Nature, 
                            levels = c("Physical",
                                       "Chemical-Physical",
                                       "Chemical",
                                       "Biological-Chemical",
                                       "Biological",
                                       "Mixed"))

light.cols <- c("#FDF2CC", "#F1F2D2", "#E2F0D9", "#DDEAE7", "#DAE3F3", "#F2F2F2")
dark.cols <- c("#FCE699", "#E1E3A7", "#C5E0B4", "#BDD3CF", "#B4C7E7", "#E4E4E4")


ggplot(data = nature_data, 
       aes(x = Year,
           y = Count, 
           group = Nature, 
           fill = Nature)) +
  geom_stream(type = "proportion", bw = 0.65, n_grid = 1000,
              color = rgb(1, 1, 1), lwd = 0.2) +
  scale_fill_manual(values = dark.cols) +
  theme_minimal() +
  theme(panel.grid = element_blank()) +
  guides(fill = FALSE) +
  labs(y = "Proportion of Stressors")

Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as of ggplot2 3.3.4.

rm(nature_data)

• Stressor classes over time

#first convert to proportions (as not using geom_stream anymore)
class_data <- class_data %>%
  group_by(Year) %>%
  mutate(proportion = Count / sum(Count))

Subset

class_data_sub <- class_data %>%
  filter(value %in% c("acidity", "habitat alteration", "metals",
                       "temperature", "UV light", "nanoparticles"))

class_data_sub$value <- factor(class_data_sub$value,
                               levels = desired_order)

ggplot(class_data_sub, 
       aes(x = Year,
           y = proportion)) +
  geom_smooth(method = "loess", span = 1,
              color = rgb(0.5, 0.5, 0.5), se = F) +
  geom_point(size = 0.5, alpha = 0.3) +
  theme_minimal() +
  labs(x = NULL, y = NULL) +
  theme(panel.background = element_rect(fill = "white", colour = "grey50")) +
  theme(panel.grid = element_blank()) +
  theme(axis.text.y = element_text(size = 4)) +
  facet_wrap(~value, scales = "free_y")  

NA
NA

All others (with 20+ occurrences)

class_data <- class_data %>%
  filter(!value %in% c("sound", "radiation", "domestic species", 
                       "biological control",
                       
                       "nanoparticles", "temperature", "UV light", 
                       "acidity", "habitat alteration", "metals"))

class_data$value <- factor(class_data$value, levels = desired_order)

ggplot(class_data, 
       aes(x = Year,
           y = proportion)) +
  geom_smooth(method = "loess", span = 1, 
              color = rgb(0.75, 0.75, 0.75), se = F) +
  geom_point(size = 0.5, alpha = 0.3) +
  theme_minimal() +
  labs(x = NULL, y = NULL) +
  theme(panel.background = element_rect(fill = "white", colour = "grey50")) +
  theme(panel.grid = element_blank()) +
  facet_wrap(~value, scales = "free_y", nrow =3, ncol =7)  

rm(class_data, class_data_sub)

• Co-occurrence of stressor classes

combinations_sub <- combinations 

# Select class columns
class_network <- combinations_sub %>%
  select(contains("Class"))

# Get unique classes
unique_classes <- unique(unlist(class_network, use.names = FALSE), na.rm = TRUE)

# Initialize an empty co-occurrence matrix
co_occurrence_matrix <- matrix(0, nrow = length(unique_classes), ncol = length(unique_classes),
                                dimnames = list(unique_classes, unique_classes))

# Iterate through each row of the data frame
for (i in 1:nrow(class_network)) {
  # Get the classes in the current row
  row_classes <- na.omit(unlist(class_network[i, ]))
  
  # Increment co-occurrence counts for all pairs of classes in the row
  if (length(row_classes) > 1) {
    class_pairs <- combn(row_classes, 2)
    for (j in 1:ncol(class_pairs)) {
      class1 <- class_pairs[1, j]
      class2 <- class_pairs[2, j]
      co_occurrence_matrix[class1, class2] <- co_occurrence_matrix[class1, class2] + 1
      co_occurrence_matrix[class2, class1] <- co_occurrence_matrix[class2, class1] + 1
    }
  }
}

# Divide the diagonal cells by 2 (pairs of the same class have been counted twice)
diag(co_occurrence_matrix) <- diag(co_occurrence_matrix) / 2

# create full empty matrix with my desired order (based on nature of classes in the taxonomy)
full_matrix <- matrix(0, nrow = 31, ncol = 31)
colnames(full_matrix) <- rownames(full_matrix) <- desired_order

# Loop through the row and column names of the co_occurrence_matrix
for (row_name in rownames(co_occurrence_matrix)) {
  for (col_name in colnames(co_occurrence_matrix)) {
    # Check if the row and column names exist in the full_matrix
    if (row_name %in% rownames(full_matrix) && col_name %in% colnames(full_matrix)) {
      # Add the value from co_occurrence_matrix to the corresponding cell in full_matrix
      full_matrix[row_name, col_name] <- full_matrix[row_name, col_name] + co_occurrence_matrix[row_name, col_name]
    }
  }
}


##### Some manual checks to make sure everything is working as expected ######

### check rows with at least two of a class
#filtered_data <- class_network %>%
#  filter(rowSums(. == "salinity", na.rm = TRUE) >= 2)

### check rows with at least one of a class 
#filtered_data <- class_network %>%
#  rowwise() %>%
#  filter(any(c_across(everything()) == "radiation")) %>%
#  ungroup()

rm(class_pairs, class1, class2, i, j, row_classes, 
   class_network, co_occurrence_matrix)

co_occur_class <- full_matrix

Figure 3 - Stressor Identities

## filter identities (only if they occur five times)
identity_data <- experiments %>%
  select(contains("Identity")) %>%
  gather(value = "Identity")
class_data <- experiments %>%
  select(contains("Class")) %>%
  gather(value = "Class")
taxonomy <- data.frame(Class = class_data$Class, Identity = identity_data$Identity)
taxonomy <- taxonomy %>%
  group_by(Class,Identity) %>%
  summarize(Count = n()) %>%
  filter(Count >= 5) %>%
  drop_na() %>%
  mutate(Nature = ifelse(Class %in% c("habitat alteration", "hydrology",  
                   "radiation", "sound", "temperature",
                   "UV light", "visible light", "water clarity"),
                   "Physical", 
                   ifelse(Class %in% c("biological alterations", 
                                       "biological control", "disease",
                                       "domestic species", "non-native species"),
                   "Biological",
                   ifelse(Class == "composite stressors", 
                          "Mixed", 
                        ifelse(Class == "cyanotoxin", "Biological-Chemical",
                               ifelse(Class %in% c("microplastics", "nanoparticles"), 
                                      "Chemical-Physical",
                                      "Chemical")))))) %>%
  
  mutate(Colour = ifelse(Nature == "Physical", "#FCE699",
                         ifelse(Nature == "Chemical", "#C5E0B4", 
                                ifelse(Nature == "Biological", "#B4C7E7", 
                                       ifelse(Nature == "Mixed", "#c4c4c4",
                                              ifelse(Nature == "Biological-Chemical",
                                                     "#BDD3CF", "#E1E3A7")))))) %>%
  mutate(Labels = ifelse(Count >=15, Identity, ""))

`summarise()` has grouped output by 'Class'. You can override using the `.groups` argument.

  
combinations_sub <- combinations #%>%
  #filter(Year %in% c(1965:2021))

# Select class columns
identity_network <- combinations_sub %>%
  select(contains("Identity")) %>%
  mutate(across(everything(), ~ifelse(. %in% taxonomy$Identity, ., NA)))
  

# Get unique classes
unique_identities <- unique(unlist(identity_network, use.names = FALSE), na.rm = TRUE)

# Initialize an empty co-occurrence matrix
co_occurrence_matrix <- matrix(0, nrow = length(unique_identities), 
                               ncol = length(unique_identities),
                                dimnames = list(unique_identities, 
                                                unique_identities))

# Iterate through each row of the data frame
for (i in 1:nrow(identity_network)) {
  # Get the classes in the current row
  row_identities <- na.omit(unlist(identity_network[i, ]))
  
  # Increment co-occurrence counts for all pairs of classes in the row
  if (length(row_identities) > 1) {
    identity_pairs <- combn(row_identities, 2)
    for (j in 1:ncol(identity_pairs)) {
      identity1 <- identity_pairs[1, j]
      identity2 <- identity_pairs[2, j]
      co_occurrence_matrix[identity1, identity2] <- co_occurrence_matrix[identity1, identity2] + 1
      co_occurrence_matrix[identity2, identity1] <- co_occurrence_matrix[identity2, identity1] + 1
    }
  }
}

# Divide the diagonal cells by 2 (pairs of the same class have been counted twice)
diag(co_occurrence_matrix) <- diag(co_occurrence_matrix) / 2

rm(identity_pairs, identity1, identity2, i, j, 
   row_identities, identity_network, unique_identities)

# Get rid of the row that has a name that is NA
filtered_matrix <- co_occurrence_matrix[-2, -2]